#!/bin/bash
RUNROOT=`pwd`
trajdir="/home/rhome/wangyong/zhaolingci_io4/Go-DCL1-A/kinetics"
Ndxfile="/home/rhome/wangyong/zhaolingci_io4/Go-DCL1-A/thermodynamics/DCL1-A.ndx"
Qcal="/home/rhome/wangyong/zhaolingci_io4/Go-DCL1-A/ana/thermo"
GROMACSDIR="/home/rhome/wangyong/opt/GROMETA4.5.5b/bin"
 if [ $# -lt 5 ]; then
 echo "$0 Filename TemperatureScale"
 echo "$0 1 200 1 0 0"
 echo "#mkdir  tmpdir"
 echo "$0 1 200 0 1 0"
else
    one=$1
    two=$2
    cal=$3
    wham=$4
    plot=$5
fi

tmpdir="/home/rhome/wangyong/zhaolingci_io4/Go-DCL1-A/ana/tmp/Kin$one"
 mkdir $tmpdir
for first in `seq 12 1 $one`
do
for second in `seq 1 1 $two`
do
if [ $cal -eq 1 ]; then
   cd $trajdir/${first}_${second}
cat >kinetic.SGE <<EOF
#$ -S /bin/bash
#$ -cwd
#$ -V
#$ -N K${first}-${second}
##==============Energy
echo "43 41 42 40 0" | $GROMACSDIR/g_energy -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/ener.edr  -o $tmpdir/ener$second.xvg
awk 'NR>22 {print \$2,\$3,\$4,\$5}' $tmpdir/ener$second.xvg > $tmpdir/ener$second.dat
##==============Gyrate
echo "1" | $GROMACSDIR/g_gyrate -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/traj.xtc  -o $tmpdir/gyrate$second.xvg 
awk 'NR>22 {print \$2}' $tmpdir/gyrate$second.xvg > $tmpdir/gyrate$second.dat
##==============Distance
echo "1 11" | $GROMACSDIR/g_dist -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/traj.xtc -n $Ndxfile -o $tmpdir/dist$second.xvg
awk 'NR>22 {print \$2}' $tmpdir/dist$second.xvg > $tmpdir/dist$second.dat
####===========Distance_alpha1-RNA
echo "11 13" | $GROMACSDIR/g_dist -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/traj.xtc -n $Ndxfile -o $tmpdir/dista1$second.xvg
awk 'NR>22 {print \$2}' $tmpdir/dista1$second.xvg > $tmpdir/dista1$second.dat
####===========Distance-alpha2-RNA
echo "11 14" | $GROMACSDIR/g_dist -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/traj.xtc -n $Ndxfile -o $tmpdir/dista2$second.xvg
awk 'NR>22 {print \$2}' $tmpdir/dista2$second.xvg > $tmpdir/dista2$second.dat
##==============Q
echo "0" | $GROMACSDIR/trjconv -s $trajdir/$first\_$second/$second.tpr  -f $trajdir/$first\_$second/traj.xtc -o $tmpdir/tmp$second.pdb -b 1
# $Qcal/QCAL.x $Qcal/ca/binding.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 $Qcal/QCAL.x $Qcal/binding73.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/Qplub$second.dat
 $Qcal/QCAL.x $Qcal/binding73.list 1 2 3 $tmpdir/tmp$second.pdb 196 Q 0.3
 mv Qnorm.dat $tmpdir/Qnorm$second.dat
 $Qcal/QCAL.x $Qcal/helixOne.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/QpluhO$second.dat
 $Qcal/QCAL.x $Qcal/helixSecond.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/QpluhS$second.dat
 $Qcal/QCAL.x $Qcal/sheet.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/Qplus$second.dat
## cat $tmpdir/Qplub$second.dat>> $tmpdir/Qball.dat
# $Qcal/QCAL.x $Qcal/ca/foldingPROTEIN.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 $Qcal/QCAL.x $Qcal/folding_csu.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/Qpluf$second.dat
 $Qcal/QCAL.x $Qcal/alpha1.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/QpluaO$second.dat
 $Qcal/QCAL.x $Qcal/alpha2.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/QpluaS$second.dat
 $Qcal/QCAL.x $Qcal/beta2beta3.list 1 2 3 $tmpdir/tmp$second.pdb 196 P 0.3
 mv Qplu.dat $tmpdir/QpluBeta$second.dat
## cat $tmpdir/Qpluf$second.dat>>$tmpdir/Qfall.dat
#==============Kinetics
 paste $tmpdir/Qplub$second.dat $tmpdir/Qpluf$second.dat $tmpdir/QpluhO$second.dat $tmpdir/QpluhS$second.dat $tmpdir/Qplus$second.dat $tmpdir/QpluaO$second.dat $tmpdir/QpluaS$second.dat $tmpdir/QpluBeta$second.dat $tmpdir/Qnorm$second.dat |awk '{printf("%6.3f %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n",\$1,\$2,\$4,\$6,\$8,\$10,\$12,\$14,\$16,\$18)}'> $tmpdir/k$second.dat
# $Qcal/kinetics.x $tmpdir/k$second.dat
 /home/rhome/wangyong/zhaolingci_io4/Go-DCL1-A/ana/kinetic/FPT.x $tmpdir/k$second.dat &>/dev/null

rm $tmpdir/ener$second.xvg
rm $tmpdir/gyrate$second.xvg
rm $tmpdir/dist$second.xvg
rm $tmpdir/dista1$second.xvg
rm $tmpdir/dista2$second.xvg
EOF
  qsub kinetic.SGE
fi
done
done

if [ $wham -eq 1 ]; then
for second in `seq 1 1 $two`
do
# cat $tmpdir/Qplub$second.dat>> $tmpdir/Qball.dat
# cat $tmpdir/Qnorm$second.dat>> $tmpdir/Qnormall.dat
# cat $tmpdir/Qpluf$second.dat>>$tmpdir/Qfall.dat
# cat $tmpdir/QpluhO$second.dat>>$tmpdir/QhOall.dat
# cat $tmpdir/QpluhS$second.dat>>$tmpdir/QhSall.dat
# cat $tmpdir/Qplus$second.dat>>$tmpdir/Qsall.dat
# cat $tmpdir/QpluaO$second.dat>>$tmpdir/QaOall.dat
# cat $tmpdir/QpluaS$second.dat>>$tmpdir/QaSall.dat
# cat $tmpdir/QpluBeta$second.dat>>$tmpdir/QBetaall.dat
# cat $tmpdir/ener$second.dat>>$tmpdir/Energy.dat
# cat $tmpdir/gyrate$second.dat>>$tmpdir/Gyrate.dat
# cat $tmpdir/dist$second.dat>>$tmpdir/Dist.dat
# cat $tmpdir/dista1$second.dat>>$tmpdir/Dist1.dat
# cat $tmpdir/dista2$second.dat>>$tmpdir/Dist2.dat
 echo " $second ##"|paste - $trajdir/${first}_${second}/OUTPUT >>$tmpdir/result.output
# echo "$second"|paste - $trajdir/${first}_${second}/OUTPUT100 >>$tmpdir/result.output100
done
# paste $tmpdir/Qball.dat $tmpdir/Qfall.dat $tmpdir/QhOall.dat $tmpdir/QhSall.dat $tmpdir/Qsall.dat $tmpdir/QaOall.dat $tmpdir/QaSall.dat $tmpdir/QBetaall.dat $tmpdir/Qnormall.dat|awk '{printf("%6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f  %6.3f\n",$2,$4,$6,$8,$10,$12,$14,$16,$18)}'>Qplot.dat
 cp Qplot.dat total$one.dat
 rm Qplot.dat
#cat >Qplots.SGE <<EOF
##$ -S /bin/bash
##$ -cwd
##$ -V
##$ -N Qplots${first}
# ~/opt/wham/wham2b Qplot.dat 1 2 100 100 # F(Q_b,Q_f)
# ~/opt/wham/wham2b Qplot.dat 9 2 100 100 # F(Q_b,Q_f)
# ~/opt/wham/wham2b Qplot.dat 1 3 100 100 # F(Q_b,Q_hO)
# ~/opt/wham/wham2b Qplot.dat 1 4 100 100 # F(Q_b,Q_hS)
# ~/opt/wham/wham2b Qplot.dat 1 5 100 100 # F(Q_b,Q_s)
# ~/opt/wham/wham2b Qplot.dat 3 4 100 100 # F(Q_hO,Q_hS)
######folding#######
# ~/opt/wham/wham2b Qplot.dat 2 6 100 100 # F(Q_f,Q_aO)
# ~/opt/wham/wham2b Qplot.dat 2 7 100 100 # F(Q_f,Q_aS)
# ~/opt/wham/wham2b Qplot.dat 2 8 100 100 # F(Q_f,Q_Beta)
#EOF
# qsub Qplots.SGE
#~/opt/wham/wham1 Qball.dat 2 # F(Q_Aspe)
#~/opt/wham/wham1 Qbnormall.dat 2 # F(Q_Ospe)
#~/opt/wham/wham1 Qfall.dat 2 # F(Q_Aspe)
#~/opt/wham/wham1 Qfnormall.dat 2 # F(Q_Ospe)
#~/opt/wham/wham2b Qplot.dat 1 2 100 100 # F(Q_Aspe,Q_Ospe)
#~/opt/wham/wham2b QPLU.smth 6 9 100 100 # F(Q_C0L3H,Q_atp)
fi

if [ $plot -eq 1 ]; then
cat <<EOF >plotQ.bat
#!/usr/bin/gnuplot
set encoding iso_8859_1
set term pdf enhanced color solid size 4in, 4in

#============================
#F(Q_LC), F(Q_NC), F(Q_LN)
#============================
#set xlabel "Time Steps"
#set x2label "Free Energy (K_bT)"
#set x2tics
#set ylabel "Q_{Aspe}"
#set output "FQAspe.pdf"
#p 'QPLU.smth' u 1:2 w l ti "Q_{Aspe}",\
#  'QPLU.smth.wam2' u 2:1 w l axis x2y1 lw 4 lt 7 ti "F({Q_{Aspe}})"
#============================
#F(Q_LC,Q_NC), F(Q_LC,Q_ATP), F(Q_NC,Q_AMP)
#============================
set pm3d at bs
set palette color
#set palette model CMY
set palette defined ( 0 "black", 1 "blue", 2 "cyan", 3 "green", 4 "yellow", 5 "orange", 6 "coral", 7 "red")
set view map
#set cblabel 'Energy Bar' textcolor lt 3
set cbrange [-12:0]

set title "F(Q_f,Q_i)"
set xlabel "Q_{f}"
set ylabel "Q_{i}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FreeEnergy${first}.pdf"
splot "Qplot.dat.wam12" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_f,Q_i)"
set xlabel "Q_{f}"
set ylabel "Q_{i}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FreeEnergy${first}dis.pdf"
splot "Qplot.dat.wam92" using 2:1:3 w p pt 5 ps 0.3 palette noti


set title "F(Q_{helixOne},Q_i)"
set xlabel "Q_{helixOne}"
set ylabel "Q_{i}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEhO${first}.pdf"
splot "Qplot.dat.wam13" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{helixSecond},Q_i)"
set xlabel "Q_{helixSecond}"
set ylabel "Q_{i}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEhS${first}.pdf"
splot "Qplot.dat.wam14" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{loop2},Q_i)"
set xlabel "Q_{loop2}"
set ylabel "Q_{i}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEs${first}.pdf"
splot "Qplot.dat.wam15" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{hS},Q_{hO})"
set xlabel "Q_{hS}"
set ylabel "Q_{hO}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEHOS${first}.pdf"
splot "Qplot.dat.wam34" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{alphaOne},Q_f)"
set xlabel "Q_{alphaOne}"
set ylabel "Q_{f}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEaO${first}.pdf"
splot "Qplot.dat.wam26" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{alphaSecond},Q_f)"
set xlabel "Q_{alphaSecond}"
set ylabel "Q_{f}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEaS${first}.pdf"
splot "Qplot.dat.wam27" using 2:1:3 w p pt 5 ps 0.3 palette noti

set title "F(Q_{Beta2&3},Q_f)"
set xlabel "Q_{Beta}"
set ylabel "Q_{f}"
set xrange [-0.04:1.0]
set yrange [-0.04:1.0]
set output "FEBeta${first}.pdf"
splot "Qplot.dat.wam28" using 2:1:3 w p pt 5 ps 0.3 palette noti

EOF
/home/rhome/wangyong/software/Gnuplot4.5/bin/gnuplot plotQ.bat
rm  Qplot*
fi
